با استفاده از داده های زلزله ها در ایران و جهان به سوالات زیر پاسخ دهید.

۱. با استفاده از داده های historical_web_data_26112015.rds و استفاده از نمودار پراکنش سه بعدی بسته plotly نمودار طول، عرض و عمق زلزله ها را رسم نمایید. علاوه بر آن بزرگی هر نقطه را برابر بزرگی زمین لرزه قرار دهید. 

sequake = read_rds("week_11/data/historical_web_data_26112015.rds")
plot_ly(x=sequake$Latitude,y=sequake$Longitude,z=sequake$Depth,type="scatter3d",mode='markers',size=sequake$Magnitude, color=sequake$Magnitude)

۲. پویانمایی سونامی های تاریخی را بر حسب شدت بر روی نقشه زمین رسم نمایید.(از داده زلزله های بزرگ استفاده نمایید.) 

disaster = read_delim("week_11/data/disaster.txt", "\t", escape_double = FALSE, trim_ws = TRUE)

disaster %>% filter(YEAR>2000) %>% 
rename(lat = LATITUDE,lon = LONGITUDE, name = COUNTRY,sequence = YEAR) %>%
  mutate(z=EQ_PRIMARY/100, color = colorize(FOCAL_DEPTH)) %>% 
  select(lat, lon, z, color, name, sequence) -> dis 

hcmap() %>% 
  hc_add_series(data = dis, type = "mapbubble", maxSize = "10")  %>% 
  hc_plotOptions(series = list(showInLegend = FALSE))

۳. نمودار چگالی دو بعدی زلزله های تاریخی ایران را رسم کنید.( از داده iran_earthquake.rds و لایه stat_density_2d استفاده نمایید). 

irequake = read_rds("week_11/data/iran_earthquake.rds")

#irmap = get_map("Iran",zoom = 5, color="bw")
irmap = read_rds("week_11/data/my_map.rds")

ggmap(irmap, extent = "panel", maprange=FALSE) +
   geom_density2d(data = irequake, aes(x = Long, y = Lat)) +
   stat_density2d(data = irequake, aes(x = Long, y = Lat, fill = ..level.., alpha = ..level..),
                  size = 0.01, bins = 16, geom = 'polygon') +
   scale_fill_gradient(low = "green", high = "red") +
   scale_alpha(range = c(0.00, 0.25), guide = FALSE) +
   theme(legend.position = "none", axis.title = element_blank(), text = element_text(size = 12))

۴. احتمال اینکه در ایران در پنج سال آینده زلزله به بزرگی هفت ریشتر رخ دهد را محاسبه کنید. (از احتمال شرطی استفاده کنید.)

۵. بر اساس داده های زلزله های بزرگ ابتدا تعداد و متوسط کشته زلزله ها را بر حسب کشور استخراج نمایید. سپس نمودار گرمایی تعداد کشته ها را بر روی کره زمین رسم نمایید.(مانند مثال زیر!) 

disaster %>% group_by(COUNTRY) %>% 
  summarise(mean_deaths = round(mean(TOTAL_DEATHS, na.rm = T)), all_deaths = sum(TOTAL_DEATHS, na.rm = T)) -> death_data

death_data$COUNTRY = str_c(str_sub(death_data$COUNTRY, 1,1), str_to_lower(str_sub(death_data$COUNTRY, 2,-1), locale = "en"))

mapdata = get_data_from_map(download_map_data(url = "custom/world.js"))

hcmap(url = "custom/world.js", data = death_data, value = "mean_deaths",
      joinBy = c("name", "COUNTRY"), name = "mean death",
      dataLabels = list(enabled = TRUE, format = '{point.name}'),
      borderColor = "#FAFAFA", borderWidth = 0.1)
hcmap(url = "custom/world.js", data = death_data, value = "all_deaths",
      joinBy = c("name", "COUNTRY"), name = "total death",
      dataLabels = list(enabled = TRUE, format = '{point.name}'),
      borderColor = "#FAFAFA", borderWidth = 0.1)

۶. با استفاده از داده لرزه های بزرگ و به وسیله طول، عرض، شدت، عمق مدلی برای پیش بینی تعداد کشته های زلزله بیابید. 

data = disaster %>% select(LATITUDE, LONGITUDE, FOCAL_DEPTH, EQ_PRIMARY, DEATHS)
data = na.omit(data)

sample_size = floor(0.8*nrow(data))
train_index = sample(seq_len(nrow(data)), sample_size, replace = FALSE, prob = NULL)
train_data = data[train_index,]
test_data = data[-train_index,]

trainMod = lm(DEATHS ~ LATITUDE + LONGITUDE + FOCAL_DEPTH + EQ_PRIMARY, data = train_data)
summary(trainMod)

Call:
lm(formula = DEATHS ~ LATITUDE + LONGITUDE + FOCAL_DEPTH + EQ_PRIMARY, 
    data = train_data)

Residuals:
   Min     1Q Median     3Q    Max 
-11904  -3244  -1728     78 311835 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -14201.885   3314.196  -4.285 2.01e-05 ***
LATITUDE        58.387     23.178   2.519   0.0119 *  
LONGITUDE       -2.834      6.115  -0.463   0.6432    
FOCAL_DEPTH    -19.639     11.321  -1.735   0.0831 .  
EQ_PRIMARY    2476.990    496.877   4.985 7.37e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14640 on 941 degrees of freedom
Multiple R-squared:  0.02878,   Adjusted R-squared:  0.02465 
F-statistic: 6.971 on 4 and 941 DF,  p-value: 1.568e-05
train_data$pred = predict(trainMod, train_data, type = 'response') 
head(train_data)
# A tibble: 6 x 6
  LATITUDE LONGITUDE FOCAL_DEPTH EQ_PRIMARY DEATHS   pred
     <dbl>     <dbl>       <int>      <dbl>  <int>  <dbl>
1     1.60     -79.4          33        7.7    600  4541.
2   -16.7      -72.7          87        6.8      1   166.
3    -7.96     110.           13        6.3   5749   370.
4    11.9      -86.0          33        5.4      7  -537.
5    27.4       88.7          10        3.5      2 -4377.
6    40.7       30.3          22        5.7      2  1778.

۷. با استفاده از داده worldwide.csv به چند سوال زیر پاسخ دهید. تحقیق کنید آیا می توان از پیش لرزه، زلزله اصلی را پیش بینی کرد؟

برای این سوال از مدل داده شده در دیتای استفاده شده استفاده کرده ام.

و به صورت مستقل مدلی فیت نکرده ام. 

wquake = read.csv('week_11/data/worldwide.csv')
print('number of predicted eqrchquakes / all earthquakes')
[1] "number of predicted eqrchquakes / all earthquakes"
sum(!is.na(wquake$rms)) / length(wquake$rms)
[1] 0.999901
print('mean rms is')
[1] "mean rms is"
mean(wquake$rms, na.rm = TRUE)
[1] 0.683795
print('min mag is')
[1] "min mag is"
min(wquake$mag, na.rm=TRUE)
[1] 2.5
print('max mag is')
[1] "max mag is"
max(wquake$mag, na.rm=TRUE)
[1] 8.2

۸. گزاره " آیا شدت زلزله به عمق آن بستگی دارد" را تحقیق کنید؟ (طبیعتا از آزمون فرض باید استفاده کنید.)

با توجه به پی-ولیوو شدت زلزله با عمق آن اندکی همبستگی دارد. 

cor(wquake$mag, wquake$depth)
[1] 0.1991147
cor.test(wquake$mag, wquake$depth)

    Pearson's product-moment correlation

data:  wquake$mag and wquake$depth
t = 50.03, df = 60629, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1914584 0.2067469
sample estimates:
      cor 
0.1991147

۹. میانگین سالانه زلزله ها را بر حسب کشور به دست آورید. آیا میتوان دلیلی در تایید یا رد تئوری هارپ ارائه کرد. 

wquake$place = as.character(wquake$place)

wquake %>% filter(str_detect(place, ', ')) -> new_wquake
x = unlist(lapply(strsplit(new_wquake$place, ', '), '[[', 2 ))

new_wquake %>% mutate(country = x, year=substr(time, 1, 4)) -> new_wquake

new_wquake %>% group_by(country, year) %>% summarise(earthquake_cnt=n()) -> mean_earthquakes

head(mean_earthquakes)
# A tibble: 6 x 3
# Groups:   country [2]
  country     year  earthquake_cnt
  <chr>       <chr>          <int>
1 Afghanistan 2015              15
2 Afghanistan 2016             193
3 Afghanistan 2017             159
4 Afghanistan 2018              46
5 Alabama     2016               4
6 Alabama     2017               6

۱۰. سه حقیقت جالب در مورد زلزله بیابید.

استان های زلزله خیز ایران 

# Top 6 Provinces in Iran with most earthquakes.
sequake %>% group_by(Province) %>% summarise(count = n()) %>% arrange(desc(count))-> result
head(result)
# A tibble: 6 x 2
  Province        count
  <chr>           <int>
1 kerman            169
2 Khorasan Razavi   166
3 isfahan           157
4 Semnan            141
5 Yazd              132
6 East Azarbaijan   130

کشور های زلزله خیز دنیا 

# Top 6 countries in the world with most earthquakes.
disaster %>% group_by(COUNTRY) %>% summarise(count = n()) %>% arrange(desc(count))-> result
head(result)
# A tibble: 6 x 2
  COUNTRY   count
  <chr>     <int>
1 CHINA       585
2 JAPAN       401
3 IRAN        373
4 INDONESIA   368
5 ITALY       325
6 TURKEY      321